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ABSTRACT 

We  examined  the  effects  of  subdural  electrical  stimulation  on  a  high-density  network  consisting  of  several  populations  of  multicompartment 
cell  types.  The  results  can  be  summarized  as  follows:  1)  Electrical  stimulation  mainly  affects  and  activates  axon  initial  and  the  most  distal 
axonal  segments  in  neurons.  The  most  distal  axonal  segments  are  locations  where  presynaptic  action  potentials  can  originate  regardless  of 
the  state  of  the  axon  initial  segment  (i.e.  depolarization,  hyperpolarization)  2)  The  effect  of  dendritic  arbor  structure  on  the  axonal  activation 
threshold  is  most  prominent  in  the  case  of  multipolar  neurons  with  large-diameter  basal/apical  dendrites  that  are  oriented  parallel  to  the 
electric  field  lines.  3)  The  timing  of  presynaptic  terminal  activation  in  neurons  subjected  to  electrical  stimulation  is  not  solely  determined  by 
the  axonal  delay  (i.e.  orthodromic  propagation)  but  rather  depends  on  the  details  of  the  applied  stimulation  field,  axonal  branching  structure 
and  axon  orientation  in  respect  to  the  electrode  position.  4)  A  single  stimulation  pulse  causes  a  sequence  of  action  potentials  ectopically 
generated  in  axons,  which  in  turn  produce  a  temporal  variation  in  the  timing  of  postsynaptic  neurons  activation. 
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Scientific  Progress 


(1)  Forward 


This  is  the  final  report  on  ARO  grant  W91  INF-1 2-1 -041 8  titled  "Biophysical  model  of  cortical  network  activity  and  the  influence 
of  electrical  stimulation”.  We  present  a  record  of  our  activities  during  this  three-year  project  dated  15  August  2012  -  14  August 
2015.  This  project  was  focused  on  the  developing  of  new  methods  and  modeling  techniques  allowing  for  the  simulation  of  active 
networks  subjected  to  electrical  stimulation.  The  report  describes  modeling  studies  that  were  carried  out  as  an  attempt  to  model 
and  quantify  the  effects  of  subdural  electrical  stimulation  on  a  high-density  network  consisting  of  several  populations  of 
multicompartment  cell  types. 

(2)  Table  of  Contents 
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(6) Bibliography 

(7) DD  Form  882 

(3)  Statement  of  Problem  Studied 

This  grant  aimed  to  use  a  biologically  realistic  network  model  to  investigate  how  electrical  stimulation  affects  behaviors  of  the 
brain’s  cortical  network  exhibiting  pathologic  or  physiologic  cortical  firing  patterns.  Modeling  efforts  were  focused  on  providing 
information  in  a  quantitative  manner  on  the  effects  of  stimulation  on  cortical  network  activity  and  its  temporal  evolution  in  the 
areas  adjacent  and  remote  to  the  stimulating  electrode.  The  onset  of  focal  seizure  was  modeled  in  in  the  built  network  model 
and  the  results  derived  from  these  modeling  studies  were  compared  with  recordings  obtained  from  subdural  microwire  arrays 
(for  local  filed  potentials  measurements)  in  humans  at  the  onset  of  the  focal  seizures. 

(4)  Summary  of  Work  Completed 

Aim  1)  To  develop  a  calibrated  computational  model  of  normal  and  epileptic  cortex  using  multicompartment  neuronal 
representations  and  an  applied  external  electric  field. 

This  specific  aim  has  been  completed.  We  developed  a  realistic  computational  model  of  cortical  network  in  which  the  effect  of 
an  external  electric  field  can  be  simulated  and  investigated.  The  developed  network  model  represents  an  unfolded  slab  of 
neocortical  tissue  (4.8  x  4.8  x3mm),  consists  of  a  several  different  types  of  multicompartmental  cortical  neurons  (excitatory  and 
inhibitory)  and  encompasses  multi-layer  structure  of  cortical  circuitry.  The  developed  model  is  unique  and  has  several  advances 
comparing  to  other  existing  models,  which  include  multi-neuronal  types  and  high  cell  density  (103/mm3).  The  model  allows  for 
better  quantification  of  the  effects  of  stimulation  in  terms  of  depth,  focality,  and  spatial  extent  of  the  recruited  neuronal  volume 
by  the  stimulating  electrode  as  well  as  for  the  creation  of  maps  of  direct  neural  activation.  In  particular,  the  developed  model 
allows  for  studying  how  the  effects  of  stimulation  vary  across  neuronal  populations  in  various  neuron  types.  It  allows  for 
identification  of  axonal  sites  where  presynaptic  action  potentials  are  initiated  in  each  neuronal  type  and  for  the  statistical 
analyses  of  these  recruited  presynaptic  action  potentials.  Based  on  these  analyses,  the  effect  of  changing  various  stimulation 
parameters  on  neuronal  recruitment  across  the  single  neuronal  population  and  types  can  be  assessed.  These  include  the  effect 
of  changing  the  amplitude  of  the  stimulating  current  and  varying  the  polarity  sequence  of  the  charge-balanced  pulse.  The  high 


density  of  neurons  in  our  model,  compared  to  densities  in  similar  previously  studied  models,  permits  the  quantification  of  how 
the  electrode  size  affects  the  number  of  recruited  neurons.  The  results  of  these  studies  have  been  published  in  the  October 
2015  issue  of  the  Neuromodulation:  Technology  at  the  Neural  Interface  journal[1]  and  were  presented  in  2014  Society  for 
Neuroscience  Meeting. 

Aim  2)  To  compare  the  results  obtained  in  our  simulation  studies  with  recordings  obtained  from  a  new  subdural  grid  design 
incorporating  standard  electrode  elements  and  microwires  for  local  field  potential  measurements  during  electrical  stimulation  for 
cortical  mapping  to  identify  context  dependent  stimulation  effects. 


This  specific  aim  has  been  completed.  We  recorded  local  field  potentials  from  microelectrode  arrays  (micro-iEEG)  with  a 
millimeter  scale  spatial  resolution  during  ictogenesis  of  focal  seizures  in  humans.  In  parallel,  a  focal  seizure  activity  was 
simulated  in  the  cortical  network  model.  The  cortical  model  yielded  analogous  simulated  local  filed  potentials  data  (9  or  16 
electrodes)  with  a  similar  millimeter  scale  spatial  resolution.  A  multivariate  autoregressive  (MVAR)  model  was  fitted  to  simulated 
and  recorded  multichannel  micro-iEEG  data  which  were  treated  as  one  multivariate  stochastic  process.  The  MVAR  model 
coefficients  were  used  to  calculate  spectral  measures  of  causality  in  the  Granger  sense  and  were  later  used  for  reconstruction 
of  the  seizure  propagation  pattern.  In  order  to  characterize  patterns  of  seizure  propagation  between  microdomains  (areas 
covered  by  microelectrodes)  we  analyzed  high  frequency  components  (high  gamma)  in  multichannel  data  obtained  from  MVAR 
model.  Analyses  in  high  gamma  range  allow  determined  dynamic  causal  interactions  between  microdomains  and  consequently 
reviled  pattern  of  micro  propagation  of  seizure  activity  inside  and  outside  of  the  seizure  epileptic  zone.  In  addition  we  compared 
high  frequency  micro  propagation  with  those  over  the  surrounding  macroelectrodes  in  low  and  high  frequency  bands.  The 
results  of  these  studies  have  been  published  in  the  August  2015  issue  of  the  Journal  of  Neural  Engineering^]  and  in  the 
proceedings  of  the  36th  Annual  International  Conference  of  the  IEEE  Engineering  in  Medicine  and  Biology  Society,  EMBC  2014 
[3], 


(5)  Summary  of  Most  Important  Results 
See  PDF  file  in  the  attachment 
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Technology  Transfer 


Our  cortical  network  model  developed  in  collaborations  with  the  Army  Research  Laboratory,  Human  Research  Directorate, 
Aberdeen  Proving  Ground  is  currently  being  used  to  investigate  the  relationships  between  blast-related  neuronal  injuries  and 
large-scale  brain  network  effects  (Dr.  Piotr  J.  Franaszczuk  and  Dr.  Jean  M.  Vettel  investigators).  These  efforts  strive  towards 
identifying  new  robust  biomarkers  in  ECoG/EEG  that  can  aid  in  determining  severity  of  the  injury,  monitoring  of  the  progression 
as  well  as  in  predicting  the  outcome.  Results  of  this  work  were  presented  in  October  201 5  at  the  Society  for  Neuroscience 
Meeting  in  Chicago,  IL. 


New  promising  applications  of  this  cortical  model  in  the  context  of  studying  relevant  topics  for  ARL  projects  are  being  currently 
discussed  with  ARL  investigators.  This  includes  a  new  collaborative  project  with  the  ARL  based  on  the  ARO  65459-LS  proposal 
for  the  development  of  a  computational  model  of  auditory  cortex  and  hearing  system  (Dr.  Dana  Boatman  and  Dr.  Anderson 
Pis).  This  project  in  addition  to  experimental  studies  heavily  relies  on  our  cortical  network  model  that  has  been  specially  tailored 
for  the  needs  of  auditory  cortical  system  processing. 


Summary  of  Most  Important  Results 


The  novelty  of  these  studies  is  twofold:  we  can  make  a  prediction  about  the  membrane 
polarization  at  each  compartment  of  each  individual  neuron  and  at  the  same  time  obtain  valuable 
statistical  information  on  the  number  of  recruited  neurons  along  with  the  initiation  sites  of  the  induced 
presynaptic  action  potentials  (PAPs).  This  approach  differs  from  earlier  studies  in  which  either  a  lower 
numbers  of  neurons  were  used  or  just  investigated  a  few  neurons  positioned  in  different  locations.  In 
contrast  to  these  previous  studies,  which  provided  valuable  but  more  qualitative  description,  our 
approach  allows  obtaining  a  detailed  quantitative  statistical  data  about  the  effects  of  subdural  electrical 
stimulation. 

Neuronal  population  recruitment 

Fig.l  illustrates  an  example  of  input-output  curves  for  the  axonal  activation  of  pyramidal  and 
basket  neurons  in  supragranular  layer  ll/lll  obtained  for  a  0.3  ms  long  charge-balanced  current  pulse 
sequence  (anodal-cathodal  vs.  cathodal-anodal).  In  the  case  of  pyramidal  neurons,  PAPs  were  initiated 
in  either  the  axon  initial  segments  (with  the  exception  of  a  low  current  amplitude  <  3  mA)  or  in  the  axon 
branch  terminal  compartments.  The  simultaneous  occurrence  of  PAPs  in  the  axon  initial  segment  (AIS) 
and  the  axon  terminal  compartment  or  in  both  terminal  compartments  was  observed  in  less  than  6%  of 
activated  axons  for  stimulating  currents  larger  than  3  mA  (anodal-cathodal  pulse).  In  contrast,  in  basket 
neurons,  the  initiation  of  PAPs  in  the  AISs  took  place  only  in  a  small  fraction  of  these  neurons  and  only  in 
the  case  of  cathodal-anodal  stimulation  (Fig.  ID).  In  the  case  of  anodal-cathodal  pulse  stimulation,  PAPs 
in  basket  neurons  were  solely  initiated  in  the  axon  terminal  compartments.  In  the  case  of  cathodal- 
anodal  pulse  stimulation,  a  strong  reduction  in  the  numbers  of  neurons  activated  in  axonal  terminals 
was  observed  in  both  pyramidal  and  basket  neurons  (Fig.  1C  and  D) 

Lateral  and  vertical  extent  of  recruited  volume 

The  depth  of  the  spatial  extent  of  neuron  recruitment  in  the  vertical  direction  (perpendicular  to 
the  electrode  surface)  varied  for  neuron  types.  Fig.  2  illustrates  depth  recruitment  curves 
(superimposed)  for  all  types  of  simulated  neurons  in  supra,  inter,  and  infragranular  layers  as  a  function 
of  the  applied  current  amplitude  to  the  electrode  for  a  0.3  ms  long  anodal-cathodal  (A)  and  cathodal- 
anodal  (B)  pulse  sequences.  Neuron  type  specific  differences  in  the  depth  of  recruitment  are  clearly 
visible  in  the  infragranular  layer  where  pyramidal  neurons  in  layers  V  and  VI  at  a  given  cortical  depth 
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Figure  1  Input-output  characteristics  of  axonal  activation  for  pyramidal  (diamonds)  and  basket  neurons  (squares)  in  the 
supragranular  layer  (ll/lll)  as  a  function  of  the  current  amplitude  (range  2-15  mA)  applied  to  the  electrode.  Legends  (left 
insets)  show  neuronal  axons  along  with  sites  of  axonal  activations  indicated  in  red.  Line  color  next  to  the  axon  matches  the 
corresponding  input-output  activation  curve  for  a  given  axonal  site.  Blue  line  corresponds  to  a  total  number  of  recruited 
neurons  (sum  of  input-output  activation  curves  for  all  axonal  sites).  A  single  0.3  ms  long  anodal-cathodal  (top  panels)  and 
cathodal-anodal  (bottom)  pulse  was  applied. 


need  approximately  30%  larger  current  amplitude  of  the  stimulating  pulse  sequence  than  basket 
neurons  in  order  to  be  activated. 

The  volume  of  recruited  neurons  under  the  electrode  has  a  radial  symmetry  and  similar  to  the 
depth  of  spatial  extend  distribution  (Fig.  2A  and  B),  the  radius  of  the  lateral  extent  of  the  recruited 
neurons  can  be  estimated  by  measuring  the  distance  between  the  center  of  the  electrode  and  the 
position  of  the  most  distal  recruited  neuron  in  planes  parallel  to  the  electrode  surface.  Fig.  2C  and  D 
show  radii  of  the  lateral  extent  of  neuron  recruitment  under  the  electrode  as  a  function  of  the  applied 
current  amplitude  for  all  simulated  neuron  types.  With  the  exception  of  pyramidal  neurons  in  the 
infragranular  layers,  there  is  a  convergence  of  curves  describing  the  radii  of  recruitment  volumes  as  the 
amplitude  of  the  applied  current  increases  (supra-  and  inter-granular  layers). 
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Figure  2  Depths  (A,B)  and  radii  (C,D)  of  a  spatial  extent  of  neuron  recruitment  under  the  electrode  as  functions  of  the  applied 
current.  Depths  and  radii  were  measured  in  directions  perpendicular  and  parallel  (on  a  plane)  to  the  pia  respectively.  Colors 
of  curves  correspond  to  simulated  neuron  in  supra,  internal  and  infragranular  layers  indicated  on  left  inset. 


Pulse  polarity  sequence 

In  terms  of  the  total  number  of  recruited  neurons  the,  anodal-cathodal  pulse  sequence  was 
found  to  be  more  efficient  than  the  equivalent  cathodal-anodal  pulse  (Fig.  2A,B  and  Fig.4A,B,C).  The 
anodal-cathodal  pulse  as  mentioned  earlier  in  general  recruited  more  neurons  in  the  axon  terminals 
than  the  cathodal-anodal  pulse.  This  can  be  seen  in  Figs.  2A  and  B  (green  line)  as  well  as  in  Figs.  4A,B  and 
C  (cadet  blue  colored  cloud)  where  in  the  case  of  the  anodal-cathodal  pulse  a  considerably  larger 
number  of  neurons  activated  in  the  axonal  terminals  was  observed  among  all  recruited  neurons.  The 
anodal-cathodal  pulse  sequence  was  also  found  more  effective  for  the  activation  of  neurons  in  deeper 
layers  than  the  analogous  cathodal-anodal  pulse.  This  can  be  seen  in  Fig.  2A  where  pyramidal  neurons  in 
layers  V  and  VI  (green  lines)  were  recruited  for  l0  >  10  mA  by  the  anodal-cathodal  pulse  while  these 
neurons  were  not  recruited  at  all  by  the  equivalent  cathodal-anodal  pulse  (Fig.  2B). 
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Figure  3  A  plot  summarizing  recruitment  of  neurons  in  layers  ll-V  in  in  a  plane  perpendicular  to  the  electrode  surface  and 
including  the  electrode  center  line.  The  volumes  corresponding  to  recruited  neurons  (blue  -  pyramidal,  turquoise  -  stellate, 
yellow  -  tufted  pyramidal  V,  red  -  basket)  may  overlap  (supragranuiar  layer)  but  in  general  are  not  equal.  The  volumes  of 
recruited  pyramidal  neurons  in  the  infragranular  layer  are  notably  smaller  than  the  volume  of  recruited  basket  neurons 
(red).  The  position  of  the  electrode  is  indicated  by  the  top  black  line.  The  applied  current  amplitude  was  l0  =  12.5mA  and  the 
stimulation  pulse  was  0.3  ms  (anodal-cathodal).  The  pyramidal  neurons  in  layer  VI  are  not  activated  at  this  particular  value 
of  l0  and  neurons  in  layer  I  were  not  simulated. 


Neuron-type  specific  recruitment 

Curves  describing  the  lateral  and  vertical  extent  of  the  recruited  neuron  volumes  (Fig.  2C  and  D) 
depend  on  neuron  type,  which  suggests  that  the  recruitment  volumes  will  not  be  identical.  The  plot 
summarizing  the  volumes  and  distribution  of  recruited  neurons  in  layers  ll-V  in  a  plane  vertical  to  the 
electrode  surface  including  the  electrode  center  line  is  shown  in  Fig.  3.  In  the  supragranuiar  layers  the 
volume  of  the  recruited  pyramidal  neurons  (blue)  overlaps  almost  exactly  with  the  volume  of  recruited 
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basket  neurons  (red)  which  is  consistent  with  Fig.  2C  for  l0  =  12.5  mA  current  amplitude  applied  to  the 
electrode.  The  recruitment  of  neurons  in  the  infragranular  layer  (V)  provides  an  example  of  unequal 
volumes.  Particularly  the  volumes  corresponding  to  the  recruitment  of  two  types  of  pyramidal  neurons 
(regular  pyramidal  -  blue  and  tufted  pyramidal  -  yellow)  lay  within  the  volume  of  the  recruited  basket 
neurons.  The  latter  however  extends  far  beyond  the  volume  of  recruited  pyramidal  neurons. 
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Figure  4  3D  maps  of  recruitment  (direct  axonal  activation)  of  pyramidal  neurons  in  the  supragranular  layer  by  a  single  0.3  ms 
pulse.  Stimulating  current  amplitudes  (3,  6  and  12  mA),  and  stimulus  polarization  along  with  the  position  of  the  electrode  are 
indicated  at  the  top  of  each  subset.  Colors  representing  the  corresponding  sites  of  axonal  initiation  (indicated  in  red)  are 
shown  on  the  right. 


Sites  of  PAP  initiation 

Besides  the  neuronal  class  specific  differences  in  neuron  recruitment,  the  response  of  a  neuron 
to  stimulation  varies  within  each  neuron  type  (Fig.  1).  The  factors  affecting  the  site  of  initiation  of  PAPs 
in  a  given  neuron  type  include  the  position  and  orientation  of  the  neuron  relative  to  the  electrode  and 
the  amplitude  of  the  applied  current  to  the  electrode.  Fig.  4  illustrates  3D  maps  showing  the  distribution 
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0.5  ms 


Figure  5  Sequence  of  excitation  of  compartments  in  a  model  of  a  basket  cell  neuron  resulting  from  a  single  pulse  of 
stimulation.  Presynaptic  spikes  are  generated  in  one  (panel  A,  left  top  insets  shows  branch  marked  in  red  color)  or 
simultaneously  in  both  distal  axonal  branches  (panel  B,  left  top  insets  branches  marked  in  red)  and  propagate  antidromically 
to  the  soma  and  later  to  the  dendritic  compartments.  The  second  PAP  occurs  at  the  axon  terminal  approximately  1.5  ms 
after  the  initial  spike  and  results  from  orthodromic  propagation  that  follows  somatic  depolarization  (yellow).  The  pulse 
timing  and  duration  are  indicated  at  the  bottom  by  the  blue  trace.  In  panel  A,  the  PAP  is  generated  in  the  left  distal  axonal 
branch  (turquoise)  and  propagates  antidromically  towards  the  distal  axon  (red)  and  next  to  the  right  axonal  branch  (green 
and  pink)  and  somatic/dendritic  compartments  (yellow/black).  Current  pulse  amplitude  in  6  mA;  duration  0.3  ms 
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of  sites  of  axonal  AP  initiation  obtained  for  pyramidal  neurons  in  supragranular  layers  and  subjected  to 
3,  6,  and  12  mA  charge-balanced  current  pulse  sequences  (anodal-cathodal  and  cathodal-anodal).  Axons 
that  are  closest  and  directly  under  the  electrode  are  predominantly  activated  in  AISs  (blue  core).  The 
initiation  within  axon  terminals  takes  place  solely  in  the  basket  like  area  (light  blue/green)  that  is 
adjacent  and  surrounds  the  core  (AIS  activation).  Simultaneous  initiation  of  PAPs  in  the  AIS  and 
terminal(s)  is  also  demonstrated  in  Fig.  4  (yellow/orange).  These  take  place  mostly  in  the  area  between 
the  core  and  the  surrounding  basket.  It  must  be  noted  that  the  activation  maps  do  not  illustrate 
instantaneous  activation  of  axonal  compartments.  The  points  used  to  create  the  maps  in  Fig.  4  were 
collected  over  a  4  ms  interval  after  the  onset  of  the  stimulation  pulse  and  might  be  associated  with 
latencies  up  to  3.1  ms  (in  the  case  of  a  cathodal-anodal  pulse).  These  studies  confirmed  that  the  sites  of 
direct  activation  of  neurons  are  AISs  but  also  indicate  that  more  distal  axon  segments  in  the  proximity  of 
the  AIS  (below  the  first  branchpoint)  are  likely  locations  where  PAPs  can  originate  regardless  of  the  state 
of  the  AIS  (sup-  or  sub-threshold  activation,  hyperpolarization). 

Impact  of  dendritic  morphology  on  axonal  activation 

After  the  stimulation,  the  level  of  polarization  of  a  given  axon  compartment  depends  on  the  net 
contribution  of  current  inputs  introduced  by  all  the  remaining  axonal  compartments.  In  the  case  of  the 
AIS,  in  addition  to  the  net  axonal  contribution,  the  proximal  dendritic  arbor  has  a  significant  impact  on 
the  polarization  of  this  compartment.  This  is  different  from  the  activation  of  the  myelinated  axon  in 
remote  nodes  of  Ranvier,  which  are  in  electrical  isolation  from  the  soma  or  the  AIS.  The  results  suggest 
that  despite  the  relatively  high  soma->axon  resistance  (small  axon  diameter)  the  mutual  somato¬ 
dendritic  polarization  may  contribute  to  the  activation  of  the  AIS.  Depending  on  the  orientation  of  a 
neuron  relative  to  the  electrode  position  and  stimulus  polarity,  different  compartments  along  the 
neuron  will  respond  with  either  a  positive  or  a  negative  polarization  of  a  different  magnitude.  For  the 
AIS  (a  perisomatic  axonal  compartment)  it  is  important  whether  the  net  axonal  and  somato-dendritic 
current  flows  make  an  antagonistic  or  synergetic  contribution.  This  in  turn  will  be  determined  by  the 
neuron  morphology,  orientation  relative  to  the  electrode  position,  and  the  stimulus  polarization.  For 
example  the  absence  of  the  AIS  activation  in  basket  neurons  in  Fig.  IB  by  an  anodal-cathodal  pulse  (for 
comparison  see  pyramidal  neurons  in  Fig.  1A)  can  be  explained  by  the  morphological  differences  in  the 
dendritic  arborization  of  basket  and  pyramidal  neurons.  The  pyramidal  neuron  model  has  large  diameter 
apical  dendrites  and  several  small  diameter  basal  dendrites.  The  basket  neuron  model,  contrary  to  the 
pyramidal  neuron,  has  four  symmetrical  (across  the  neuron  main  axis)  large  diameter  dendritic  branches 
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(Fig.  5,  top  right  insets).  For  this  stimulation  paradigm  and  the  electrode  configuration,  the  neuron's 
normal  axes  are  perpendicular  to  the  electrode  surface  meaning  that  the  electric  field  lines  directly 
under  the  electrode  remain  parallel  to  these  axes.  The  low  resistance  between  the  soma,  the  apical 
dendritic,  and  in  the  case  of  basket  neurons,  the  basal  dendritic  branches  creates  an  easy  pathway  for 
the  flow  of  axial  currents  in  the  direction  that  is  determined  by  the  stimulus  polarization.  During  the 
anodal  stimulation  the  somal  compartment  of  pyramidal  neurons  would  be  more  depolarized  because  a 
smaller  volume  of  axial  current  flows  into  the  basal  compartments.  In  the  basket  neuron,  the  somal 
compartment  is  hyperpolarized  rather  than  depolarized  during  the  same  anodal  pulse  because  the 
current  can  easily  flow  into  the  large  diameter  and  basal ly  oriented  compartment.  The  hyperpolarization 
of  the  soma  in  basket  neurons  affects  the  polarization  of  the  AIS  by  increasing  the  threshold  for 
activation  of  this  compartment.  This  may  explain  the  lack  of  PAP  in  the  AISs  in  basket  neurons  subjected 
to  an  anodal-cathodal  pulse  in  Fig  IB. 

Sealed-end  effect  in  distal  axonal  compartments 

In  the  case  of  the  most  distal  axon  compartments  (below  the  first  branchpoint),  in  addition  to 
the  net  current  contribution  from  the  adjacent  axon  compartments,  the  sealed-end  effect  will  also 
contribute  to  the  polarization  of  these  compartments.  In  this  stimulation  paradigm,  the  neuron's  normal 
axes  are  perpendicular  to  the  electrode  surface.  The  axonal  terminals  will  always  represent  more  distal 
compartments  for  the  neurons  (relative  to  the  electrode  position)  and  therefore  they  will  be  the  most 
depolarized  in  the  case  of  anodal  stimulation  and  the  most  hyperpolarized  compartments  after  cathodal 
stimulation  respectively.  This  effect  is  clearly  evident  in  Fig.  1C  and  D  where  the  number  of  neurons 
activated  in  axon  terminals  in  both  types  of  neurons  by  cathodal-anodal  pulses  is  lower  compared  to 
anodal-cathodal  pulses  (Fig.  1A,B).  There  is  no  rebound  excitation  in  these  compartments  after  the 
cathodal  phase  of  the  stimulus  because  the  width  of  the  pulse  phase  (0.15  ms)  is  too  short  to  induce 
anodal  break  excitation.  This  is  also  evident  in  Fig.  4D-F  where  the  volume  representing  activation  of 
axonal  terminals  is  notably  smaller. 

Temporal  response  of  neurons  to  stimulation 

Analyses  of  the  temporal  response  of  a  single  neuron  to  the  applied  stimulation  pulse  are 
illustrated  in  Fig.  5A  and  B.  The  temporal  sequences  of  excitation  occurring  in  neurons  are  shown  by 
plotting  voltage  traces  in  all  axonal,  somatic,  and  most  proximal  dendritic  compartments  in  two  basket 
neurons  (A  and  B)  after  a  0.3  ms  long  current  pulse.  PAPs  are  initiated  in  one  (A)  or  simultaneously  in 
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both  axonal  branch  terminals  (B).  A  single  stimulation  pulse  sequence  typically  induces  a  PAP  duplet  at 
axon  terminals,  with  the  second  PAP  arising  from  AIS  activation  by  antidromic  propagation  (towards  the 
soma)  from  the  initial  PAP.  In  panel  B,  these  duplets  (the  initial  and  the  second  PAP)  arrive  at  both 
terminals  at  the  same  time,  while  in  panel  A  there  is  latency  (0.5  ms)  in  the  arrival  of  these  spike  duplets 
at  axonal  branch  terminals.  These  results  suggests  that  after  stimulation,  the  time  of  activation  of 
postsynaptic  neurons  at  presynaptic  terminals  is  not  determined  solely  by  the  axonal  delay  (orthodromic 
propagation)  but  depends  on  details  of  the  applied  stimulation  field.  The  fact  that  after  stimulation, 

PAPs  arrive  at  the  presynaptic  terminals  of  various  branches  at  different  times  might  be  important  in 
characterizing  the  effects  and  evaluating  the  efficacy  of  electrical  stimulation.  The  temporal  variability  of 
PAPs  can  contribute  to  desynchronizaton  of  the  neural  activity  and  therefore  may  underlay 
desynchronizing  or  anti-epileptic  effects  of  stimulation. 


Electrode  diameter 

In  the  current-controlled  configuration  (constant  current)  for  the  simulations  used  in  this  work, 
the  electrode  diameter  has  a  minor  effect  on  the  recruitment  of  neurons.  Fig.  6  illustrates  input-output 
curves  of  neuronal  activation  obtained  for  three  different  electrode  diameters.  The  slope  of  the  input- 
output  curves  as  well  as  the  number  of  recruited  neurons  are  relatively  unaffected  by  the  electrode 
diameter  for  the  4  -6  mA  range  of  the  stimulation  current.  However  the  input-output  curves  diverged  for 
stimulation  currents  below  4  mA  and  when  the  applied  current  to  the  electrode  is  above  6  mA.  These 
results  suggest  that  when  the  electrode  diameter  is  changed,  the  total  number  of  recruited  neurons 
does  not  change  notably  as  long  as  the  current  l0 applied  to  the  electrode  remains  the  same  (which  is 
the  case  in  current-controlled  stimulation).  The  modest  impact  of  the  electrode  diameter  on  the  number 
of  recruited  neuron  can  be  explained  by  the  fact  that  the  lateral  and  vertical  extent  of  the  recruited 
neuron  volume  is  inversely  correlated  to  the  electrode  diameter.  In  general,  electrodes  with  larger 
diameters,  besides  providing  a  larger  lateral  extent  of  field  effects,  will  offer  a  smaller  extent  of  the 
recruitment  volume  in  the  vertical  direction.  Electrodes  with  smaller  diameters  will  offer  a  larger  vertical 
field  effect  but  with  a  smaller  lateral  extent.  The  latter  suggests  that  for  a  given  l0  smaller  electrodes  will 
generally  provide  better  recruitment  of  neurons  in  the  vertical  direction.  It  should  be  noted  that  this  is 
only  a  theoretical  supposition  because  for  a  given  amplitude  of  l0  (constant-current  stimulation),  when 
the  electrode  diameter  d  decreases,  the  V0  potential  at  the  electrode  will  increase  as  1/r  (r=d/2  is  an 
electrode  radius).  Therefore  small  diameter  electrodes  will  require  larger  V0  potentials  in  order  to 
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activate  the  same  number  of  neurons.  For  that  reason  the  simulations  described  in  this  work  with  small 
diameter  electrodes  in  situ  may  be  associated  with  a  higher  risk  of  tissue  damage  since  a  larger  current 
density  will  be  delivered  from  a  small  electrode  area. 


Pyramidal  supragranular  layer  I  I/I  II 


Figure  6  Effect  of  changing  electrode  diameter  on  numbers  of  activated  neurons.  Input-output  curves  of  axonal  activation 
for  pyramidal  neurons  in  the  supragranular  layer  (ll/lll)  vs.  the  current  amplitude  applied  to  the  electrode  (range  2.5-12.5 
mA)  are  shown  for  0.7, 1.5,  and  2.1  mm  electrode  diameters.  A  single  0.3  ms  long  anodal-cathodal  (A)  and  cathodal-anodal 
(B)  pulse  was  applied. 


Simulation  of  micro-iEEG  recording 

The  cortical  model  used  to  simulate  micro-iEEG  corresponds  to  a  square  area  of  0.8  x  0.8  mm 
sampled  by  five  or  sixteen  microelectrodes  with  the  distance  between  the  microelectrodes  0.4  mm.  A 
simulated  region  of  cortex  was  represented  by  layer  2/3  pyramidal  cells,  spaced  at  25  pm  in  each 
horizontal  direction.  The  LFP  for  the  model  is  derived  from  a  weighted  average  of  the  current  sources 
summed  over  the  cellular  compartments.  This  includes  synaptic  currents  (AMPA.  NMDA),  channel 
currents,  and  compartment  currents.  The  simulated  dataset  of  micro-iEEG  measured  from  five  virtual 
micro-electrodes  was  fitted  with  a  multivariate  autoregressive  model  (MVAR).  The  direct  transfer 
functions  (DTF),  direct  DTF  (dDTF),  partial  direct  coherence  (PDC)  and  generalized  PDC  (gPDC  )  were 
calculated  and  the  propagation  pattern  was  determined  over  2-30  Hz  range  as  shown  in  Fig  7.  The 
actual  propagation  of  the  wave  across  the  sampled  region  from  the  upper  right  corner  to  the  lower  left 
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corner  is  shown  on  the  right  side  of  the  Fig  7.  It  can  be  seen  that  although  all  four  measures  produce 
similar  flow  patterns  (initiating  at  right  upper  corner),  DTF  has  the  most  accurate  and  robust  estimation. 


Figure  7.  Determination  of  known  propagation  (right  top  to  bottom)  in  simulated  five  channel  micro-iEEG  data  using  DTF, 
dDTF,  PDC  and  gPDC  in  the  2-30  Hz  frequency  band  using  three  values  of  threshold  Ta  a  =  80,85,90. 


Micro-iEEG  recording  and  analyses 

The  micro-iEEG  recorded  from  4x4  micro-electrode  arrays  which  were  placed  within  a 
centimeter  of  the  clinical  seizure  foci  in  4  patients  with  intractable  epilepsy  were  used  to  determine  the 
seizure  propagation  pattern.  Each  patient  had  at  least  2  micro-electrode  arrays,  Ml  and  M2,  Ml 
denoting  the  one  closer  to  the  clinically  determined  seizure  focus  than  the  other  (M2)  for  most  patients. 
We  used  the  DTF  measure  to  calculate  propagation  over  the  frequency  band  determined  for  each 
patient  by  summing  over  the  DTF  values  in  that  frequency  range. 

Positions  of  the  macro-electrode  and  microelectrode  grids  are  shown  in  Fig  8.  Patients  PI  and 
P2  had  left  temporal  grids  (LFT)  with  two  micro-arrays  Ml  and  M2  as  shown  in  Fig  8.  Directions  of 
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seizure  propagation  pattern  were  determined  for  each  subject  (P1-P4).  Examples  of  these  analyses  are 
presented  below  for  subject  PI  and  they  can  be  summarized  as  follows. 


Figure  8.  The  reconstructed  brain  maps  of  the  four  patients  considered  showing  approximate  positions  of  the  macro¬ 
electrode  and  microelectrode  grids  along  with  the  boundaries  of  the  resected  region  in  red.  The  yellow  asterisks  in  PI,  P2,  P4 
indicate  the  seizure  onset  zone  macro-electrode  contacts  on  the  grid.  The  yellow  and  white  asterisks  in  P3  denote  SOZ 
macro-electrodes  during  seizures  1  and  2  respectively. 

In  PI,  the  earliest  seizure  onsets  on  LFT  were  localized  to  contacts  19  and  26  as  marked  by  yellow 
asterisks  and  the  resection  boundary  with  a  red  line  as  in  Fig.  8.  All  Ml,  M2  and  LFT  19,  22,  43,  46 
channels  were  used  to  reconstruct  propagation  across  Ml,  M2  and  the  four  macro-electrodes 
surrounding  Ml,  M2  as  shown  in  Fig.  9(A).  Fig.  9(D)  shows  the  top  four  percentile  of  the  thresholded 
DTF  values  plotted  as  propagation  across  Ml  (left  corner)  and  M2  (right  corner)  along  with  the  strongest 
flows  across  the  neighboring  macro-electrodes  in  the  4-50  Hz  range  (green)  and  70-110  Hz  range  (blue). 
The  macro  flow  over  each  time  window  was  aggregated  and  is  plotted  in  Fig.  9(B)  with  the  thickness  of 
the  arrows  being  proportional  to  the  number  of  flows  across  a  pair  of  contacts.  Both  the  frequency 
bands  have  similar  dominant  flows  except  for  the  top  contact  pair.  Furthermore,  the 
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Figure  9.  (A)  15  channels  of  Ml  (blue)  and  12  channels  of  M2  (black)  recorded  from  PI  during  a  clinical  seizure,  the  inset 
shows  a  schematic  of  the  implanted  grid  along  with  Ml  and  M2  with  the  black  asterisk  denoting  the  location  of  the  SOZ 
macro-electrode  contact  with  respect  to  Ml;  (B)  the  4-50  Hz  (green)  and  the  70-110  Hz  (blue)  flows  in  the  macro-channels 
surrounding  Ml  and  M2.  The  thickness  of  the  arrows  is  proportional  to  the  number  of  times  a  particular  direction  of  flow 
occurred  in  the  entire  interval  considered  as  in  (D);  (C)  time  evolution  and  histogram  of  micro-electrode  contact  emanating 
the  maximum  outflow  within  grids  Ml  and  M2  respectively.  The  channel  having  maximum  outflow  for  the  highest  number  of 
times  over  the  entire  seizure  duration  is  marked  in  black  and  red  for  Ml  and  M2  respectively  in  the  schematic  in  (A);  (D)  top 
four  percentile  of  propagation  in  the  70-90  Hz  frequency  band  across  Ml  (bottom)  and  M2  (top)  along  with  the  strongest 
flow  across  the  neighboring  macro-electrode  contacts  in  the  4-50  Hz  (green)  and  the  70-110  Hz  (blue)  frequency  bands.  Each 
square  element  corresponds  to  a  time  interval  of  2  s  and  consecutive  squares  are  separated  by  0.5  s.  The  temporal  evolution 
is  from  left  to  right.  The  32  patterns  shown  thus  correspond  to  a  total  duration  of  16  s  as  shown  in  (A). 
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channel  with  the  maximum  outflow  was  calculated  for  Ml  and  M2  and  is  plotted  over  time  (left)  and  as 
a  histogram  (right)  in  Fig.  9(C)  which  shows  that  in  Ml,  channel  8  acts  as  the  source  for  maximum 
number  of  time  windows  whereas  in  M2,  channels  10  and  4  have  the  highest  outflows.  From  Fig.  4(D), 
there  seems  to  be  a  more  stable  one  directional  flow  pattern  in  M2  although  a  single  electrode  (channel 
8)  dominates  outflow  in  Ml. 

Implications  for  realistic  computational  models 

We  used  a  realistic  computational  model  to  generate  microelectrode  data  with  a  known 
propagation  pattern  and  the  Granger  causality  measures  (DTF,  dDTF,  PDC,  gPDC)  to  reconstruct  the 
propagation  pattern  from  the  simulated  data.  One  of  the  major  challenges  of  building  a  realistic 
computational  model  for  LFPs  is  calibration  and  comparison  with  a  measured  signal.  Although  we  used  a 
computational  model  here  to  validate  our  propagation  reconstruction  method,  this  could  also  provide  a 
modeling  framework  for  tuning  such  computational  models  to  match  real  data  in  terms  of  propagation. 
Efforts  focused  on  reconstruction  of  pattern  of  activity  propagation  between  micro-electrodes  will 
potentially  play  an  important  role  in  constructing  and  validating  models  of  cortical  activity,  with  a 
reliance  on  human  recordings  for  calibration.  Specifically,  these  reconstructed  patterns  from  micro¬ 
electrodes  will  help  quantify  activity  spread  in  lateral  direction  similar  to  current  density  analyses,  which 
allow  for  studying  the  laminar  pattern  of  neural  activity  in  neocortex.  Therefore  details  of  reconstructed 
propagation  pattern  can  help  in  estimation  of  some  values  of  network  parameters  which  are  difficult  to 
determine  in  situ.  These  might  include  selection  of  weights,  numbers,  and  the  range  of  synaptic  contacts 
in  the  process  of  network  model  parameter  tuning  in  order  to  reproduce  the  propagation  pattern 
consistent  with  measurements  of  the  propagation  of  epileptiform  activity  in  neocortex. 
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